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We propose a computationally-efficient approach to multiscale simulation of polycrystalline ma- 
terials, based on the phase field crystal (PFC) model. The order parameter describing the density 
profile at the nanoscale is reconstructed from its slowly-varying amplitude and phase, which satisfy 
rotationally-covariant equations derivable from the renormalization group. We validate the approach 
using the example of two-dimensional grain nucleation and growth. 
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Why is it so hard to predict the properties of real mate- 
rials? Unlike simple crystalline solids, real materials, pro- 
duced by a wide range of processing conditions, contain 
defects and multiple grains that strongly impact mechan- 
ical, thermal, electrical response, and give rise to such im- 
portant phenomena as plasticity, hysteresis, work hard- 
ening and glassy relaxation. Moreover, it is frequently 
the case that a faithful description of materials process- 
ing requires simultaneous treatment of dynamics at scales 
ranging from the nanoscale up to the macroscopic. For 
example, dendritic growth, the generic mode of solidifi- 
cation of most metals and alloys, involves the capillary 
length at the nanoscale, the emergent pattern dimensions 
on the scale of microns, the thermal or particle diffusion 
length on the scale of 10~^m, in addition to the grain 
and sample size. 

Despite these obstacles, progress in rational material 
design requires a fundamental understanding of the way 
in which useful properties emerge as the mesoscale is ap- 
proached. Questions that must be addressed include: 
What is the collective behavior of assemblies of nanoscale 
objects? How best to achieve target mesoscale properties 
from nanoscale constituents? And how can the properties 
at nano-, meso- and intermediate scales simultaneously 
be captured quantitatively and predictively? 

A number of computational approaches to handle the 
range of length scales have been proposed recently jl|, Q , 
including quasi-continuum methods 3, y, l^, |Q] , the het- 
erogeneous multiscale method |2i, ^J, multi-scale molec- 
ular dynamics 0, 0, 0, multigrid variants and 
extensions of the phase field model 14]. These techniques 
strive to provide a unified description of the many scales 
being resolved, but in some cases require non-systematic 
ways to link the disparate scales to enable treatment of 
sufficiently large mesoscale systems. This can introduce 
spurious modes and excitations, and difficulties associ- 
ated with the transition between scales 1^. Most of 
this work is limited to crystalline materials with a few 
isolated defects jl5j. 

In this Letter, we propose a novel theoretical approach 



to these difficulties, by combining the phase field crys- 
tal (PFC) formalism [iM with renormalization group 
(RG)*!^, 1^ and related methods (see, e.g. |20i]), de- 
veloped for the analysis of hydrodynamic instabilities 
in spati ally -extended dynamical systems [l^, HE 1^ 
EE HE m 113 • We present effective equations at the 
mesoscale, from which the atomic density can readily be 
reconstructed, and show that this approach is capable of 
generating high fidelity representations of materials pro- 
cessing dynamics. Moreover, we demonstrate that the 
mesoscale equations — analogues of rotationally-covariant 
amplitude and phase equations in fluid convection [S^l — 
are computationally tractable and amenable in future 
work to adaptive grid techniques. 

Our approach is based on a form of the RG which uni- 
fies singular perturbation theory 2 Ij, and is a fully sys- 
tematic way to extract universal or large scale structures 
from spatially-extended dynamical systems. The basic 
idea is to start, not with a molecular dynamics model at 
the nanoscale, but with a density functional description 
(in this context, the phase field crystal model), whose 
equilibrium solutions are periodic density modulations. 
A system that is periodic at the nanoscale can be pa- 
rameterized in terms of a uniform phase and an ampli- 
tude: the amplitude describes the maximum variations 
in the density of the system through the unit cell, while 
the phase describes uniform spatial translations. A sys- 
tem with underlying periodicity, but which also contains 
defects or other nanostructure, can be represented by a 
density wave whose amplitude is at most slowly varying 
on the nanoscale, and a phase that is essentially uni- 
form everywhere, except near a defect. This observation 
suggests that the phase of the density is the appropri- 
ate dynamical variable to use for describing spatially- 
modulated nanoscale structure in a mesoscopic system, 
and in the vicinity of a defect it must be supplemented 
by the amplitude. 

For any spatially-extended pattern- forming dynamics, 
RG provides a prescription for obtaining slowly-varying 
amplitude and phase equations valid on scales much 
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larger than the nanoscale}21ll2al2al2j,l25|. These equa- 
tions possess the key advantage that their solutions are 
essentially uniform, with localized rapid variations near 
defects. The renormalization procedure used here is more 
general than real-space renormalization, which has been 
attempted in related contexts^, E^. In particular, our 
technique directly focuses on the instabilities that char- 
acterize the dynamics, which may not have simple real- 
space interpretations. Once the amplitude and phase are 
determined, the actual structure at the nanoscale (and 
above) can be reconstructed. Because the amplitude 
and phase equations describe solutions that are slowly- 
varying everywhere, except near a defect, adaptive mesh 
refinement can be used to solve the amplitude and phase 
equations. 

Phase field crystal model:- The phase field crystal 
model is a continuum, nonlinear partial differential equa- 
tion for the density p{x) that recently has been shown [iq. 
0-7] to capture realistic aspects of materials dynamics, in- 
cluding grain growth, ductile fracture, epitaxial growth, 
solidification processes, and reconstructive phase transi- 
tions. In addition, the model natively supports elasticity 
theory, both linear and non-linear, without any ad hoc 
modeling. Thus, the PFC can address the important 
problem of nanoscale strain effects, and their coupling 
from the nanoscale to the continuum. 

Let F{p} denote the coarse-grained free energy func- 
tional whose minima correspond to the equilibrium (lat- 
tice) state of a d-dimensional system, and whose corre- 
sponding chemical potential gradient drives the dynamics 
of p. A simple form of F that gives rise to a triangular 
lattice equilibrium state is the Brazovskii form 30] : 

F{p(x)} = I d^x [p (aAT + \{ql + V^)') p/2 + up^/A 

(1)" 

where a. A, Oo and u can be related to material 
properties [la. \u\ . and AT denotes the temperature dif- 
ference from some (higher) reference temperature. It 
is convenient to rewrite this free energy in dimension- 
less units, i.e., x = fqo, = p\/ul\q^^ r = aAT/Ag^, 
r = FAg^t, where F is a phenomenological constant 17], 
and F Fu/ \^q^~^ so that the equation of continuity 
for the density becomes 

d^/dt = V' ( [r + (1 + V')'] + V^') + C. (2) 

The conserved Gaussian noise will not generally be im- 
portant for describing phase transition kinetics, and so 
will henceforth be neglected here. 

The mean field phase diagram of the PFC equation 
((2j can be calculated analytically in a one mode 
approximation that is valid in the limit of small r, and 
represented in the plane of dimensionless temperature, 
r, and dimensionless average density, Three possible 
equilibrium solutions exist; a 'liquid', ifjc = di two- 
dimensional 'crystal' with triangular symmetry, tpT = 

At (^cos(gT^) cos{qTy/V^) - cos(2(7t^/\/(3))/2^ + ^/S, 
and a smectic phase which will be ignored for present 



purposes. The triangular lattice can exhibit persistent 
defect structures during the relaxation to equilibrium. 

Key to differentiating our approach from others is the 
fact that the PFC method was designed to investigate 
phenomena on diffusive time scales which are typically 
many orders of magnitude longer than the time scales 
accessible in molecular dynamics simulations. On the 
other hand, the PFC method suffers from some of the 
same drawbacks as molecular dynamics simulation: by 
resolving the nanoscale, brute force computation on a 
massive scale is required to capture mesoscale phenom- 
ena. We finesse this difficulty here by working with the 
slowly- varying amplitude and phase description. 

Mesoscale representation of the PFC:- The dynamics 
of the slowly- varying amplitude and phase describes fluc- 
tuations about a given set of lattice vectors, but must be 
covariant with respect to rotations of those lattice vec- 
tors, in order to properly describe polycrystalline ma- 
terials with arbitrarily oriented grains. A similar sit- 
uation arises in describing amplitude and phase varia- 
tions of convection rolls, and in the context of the model 
Swift-Hohenberg [sj equations, the form of the govern- 
ing equations was originally proposed by Gunaratne et al. 
[3^, and derived systematically from the RG formalism 
of Chen et al.|2ll by Grahamlla (see also ref. Q). 

The triangular phase solution is represented as 



V^(x) = ^ Aj{t) ex.p{ikj ■ x) -\- ijj^ 



(3) 



where ki = ko{—i^/3/2 — j/2), /c2 = /coj and ks = 
ko{iV^/2 — j/2) are the reciprocal lattice vectors, ko is 
the wavenumber of the pattern, i and j are unit vectors in 
the X- and ^/-directions, and Aj {j = 1,2,3) are the com- 
plex amplitude functions. The detailed derivation of the 
evolution equations for the RG form of the PFC equa- 
tion, similar to the approaches referenced above, will be 
given elsewhere. We simply present the result here due 
to space limitations 



dAi 
dt 



= C,A - 3Ai + 2IA2I' + 2IA3I') - 6V;A*A* 

(4) 

(together with appropriate permutations for A2 and A3), 
where 
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2ik, • V 



+ 2ik, • V 



(5) 

is the manifestly rotationally covariant operator. After 
solving the RG equations in Eq. the density is re- 
constructed using Eq. Q. 

Model Validation:- For ease of numerical implementa- 
tion, we have chosen to solve the amplitude equations 
about a globally flxed basis of lattice vectors. As a re- 
sult, the grain orientation information must be borne by 
the complex amplitude functions Aj. We specify this 
information through an initial condition Aj{0)^ where 
is the rotation angle measured with respect to the basis 



vectors, and the function Aj{0) is chosen such that, when 
the original field tjj is reconstructed as per Eq. the 
resulting grain is rotated by an angle 6. 

Fig. (^) shows the time evolution for the nucleation and 
growth of a two-dimensional film as predicted by the RG 
equation, starting from an initial condition of randomly- 
oriented seeds, with i/j = 0.285 and r = —0.25. The 
initial crystallite domains grow, colliding to form a poly- 
crystalline microstructure. The solutions obtained using 
the PFC equation are essentially indistinguishable from 
Fig. indicating excellent qualitative agreement. The 
key feature of the RG equation is its ability to correctly 
capture defect formation and motion. 




(a) t = 56 (b) t = 184 (c) t = 720 

FIG. 1: RG-reconstructed density at indicated times for het- 
erogeneous nucleation and growth in a 2-D film. 

As a more rigorous demonstration of accuracy. Fig. Q 
compares the grain boundary energy, 7, as a function of 
the misorientation angle 6>, predicted using the two algo- 
rithms. The initial condition for this test comprised two 
misaligned crystals separated by a narrow strip of liquid, 
on a periodic domain (see 17] for details). Fig. ((2j also 
shows the Read-Shockley equation [2^, a well known an- 
alytical result for small angle grain boundaries, that has 
been scaled to fit large misorientation data. The agree- 
ment is particularly good for low angle grain boundaries, 
and the values predicted by the RG equations closely fol- 
low the trends predicted by the PFC (from ^3) ^^e 
Read-Shockley equation. The maximum difference in the 
free energy as computed by the RG and PFC equations 
is about 1.6%. 

Computational efficiency:- Fig. 01 (inset) shows grid 
convergence behavior in the Read-Shockley test of the 
solutions to the PFC and RG equations. The crystals 
are misoriented by the maximum possible angle, 7r/6. We 
define the error eo = |||^^^||2 — ll^^lb], where ||^^^||2 is 
the 1/2 norm of the solution for a mesh spacing of Ax, 
and 11^^112 is the L2 norm obtained by Richardson ex- 
trapolation to Ax = consistent with a second order 
finite difference method. For a comparable level of ac- 
curacy, we see that Axrg ^ '^Axpfc- In a Forward 
Euler time step scheme, this leads to a stability condi- 
tion AtRG ^ 6 At PFC- Clearly, the RG equations offer 
significant opportunities for improved computational ef- 
ficiency. 

Fig. |31 compares the CPU time as a function of domain 
size Lx for the Read-Shockley test with = 3.88°, 
r = -0.25 and = 0.28, showing that the CPU time 
required for the RG equations ranges from about 5 to 
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FIG. 2: Comparison of grain boundary energy predicted by 
the RG and PFC equations, with the Read-Shockley equation. 
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FIG. 3: Scaling of CPU time versus domain length Lx for the 
PFC and RG equations. Inset shows error in the respective 
solutions with diminishing mesh spacing Ax. 



6 times less than for the PFC equations. Consistent 
with the grid convergence behavior described above, we 
chose Ax PFC = 7r/4, Atppc = 0.008, Axrg = 7r/2 and 
AtRG = 0.05. The difference in the free energy predicted 
by the RG equations and PFC equations was < 0.1%. 
Even more significantly, however, the amplitude func- 
tions can be solved on a non-uniform computational grid. 
For this problem, since the location of the boundary is 
known a priori^ it is easy to construct an appropriate non- 
uniform mesh. We chose constant Ay = 7r/2, and allowed 
Ax to vary from a minimum of 7r/2 to 27r. This reduced 
the size of the computational mesh from 257 x 257 to 
97x257. Even with this relatively naive implementation, 
we find that the speedup of the RG form compared to 
the original PFC form is close to a factor of ten (see 
Fig. El) 7 while the error in the free energy is still < 0.1%. 
We anticipate the computational benefits of solving these 
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equations on a fully adaptive mesh to be much higher. 




(a) Af (b) ^ 

FIG. 4: (a) Real component of the complex amplitude; (b) 
Density field ijj reconstructed using Eq. 0- Clockwise from 
the lower left, 6* = 0, 7r/24 and 7r/6. 

One limitation to the approach described here is that 
the orientation information causes a spatial variation in 
the amplitude when it is represented using the basis 
vectors of the triangular lattice. This is illustrated in 
Fig. (Q, showing the real part of one of the complex am- 
plitude functions, and the corresponding reconstructed 
density variable ip for t hree d ifferent orientations. The 
"beats" evident in Figure [4(a)l which contain grain orien- 
tation information, persist as the grains evolve. This phe- 



nomenon could limit the effectiveness of adaptive mesh 
methods, since the mesh has to resolve these structures. 
We note however, that if we write the complex ampli- 
tudes as Aj = ^jexp(z$j), where is the amplitude 
modulus and is the phase angle, we can formulate 
equations of motion for and V<^>j from Eq. (Q, 
fields which are uniform everywhere (no "beats" ) except 
near defects and interfaces. The resulting adaptive grid 
algorithm can thus be made to scale much more opti- 
mally Is^, with interface/grain boundary length rather 
than the area of significantly misoriented grains. We will 
present this work in a future article. 

In summary, we have shown that multiscale model- 
ing of complex polycrystalline materials microstructure 
is possible using a combination of continuum modeling 
at the nanoscale using the PFC model, RG and related 
techniques from spatially-extended dynamical systems 
theory. Our equations are well-suited for efficient adap- 
tive mesh refinement, thus enabling realistic modeling of 
large-scale materials processing and behavior. 
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